Taller de R: Estadística y programación

class: center, middle, inverse, title-slide

# Taller de R: Estadística y programación
## Lectura 13: Datos raster
### Eduard Martínez
### Universidad de los Andes | <a href="https://github.com/taller-R">ECON-1302</a>

---

name: contenido
# Ya vimos..

- [Operaciones geométricas con objetos sf](https://github.com/taller-R/clase-12)

# Hoy veremos

1. [Introducción a datos raster](#intro)

2. [Importar datos raster](#read)

3. [Operaciones geométricas con raster](#ope)

4. [Raster de varios layers](#layers)

---
class: inverse, center, middle
name: intro
# [1. Introducción a datos raster](#contenido)

---
# 1.1. Qué es un raster?

- Un raster es una matriz de celdas (o píxeles) organizadas en filas y columnas (o una cuadrícula) en la que cada celda contiene un valor que representa información.

<div align="center">
<img src="pics/raster.png" height=400>
</div>
Tomado de: [https://www.neonscience.org](https://www.neonscience.org/resources/learning-hub/tutorials/dc-raster-data-r)

---
# 1.2. Elementos de un raster

* Dimensión o numero de filas y columnas, así como el numero de bandas que tiene el raster

* La resolución de cada pixel

* Tipo de archivo

* CRS

* Información dentro de los pixelex

# 1.3. Extenciones de un raster

* .tif

* .nc

* .nc4

---
# 1.4. Algunas fuentes de información de datos raster

- **Luces nocturnas mensuales 2012-2019**

Pueden descargarse en [https://ngdc.noaa.gov](https://ngdc.noaa.gov/eog/viirs/download_dnb_composites.html) o en [https://payneinstitute.mines.edu](https://payneinstitute.mines.edu/eog/)

- **Deforestación (SIAC)**

Pueden descargarse en [http://www.siac.gov.co](http://www.siac.gov.co)

- **Deforestación (Hansen)**

Pueden descargarse en [https://data.globalforestwatch.org](https://data.globalforestwatch.org/datasets/14228e6347c44f5691572169e9e107ad)

- **GES-DISC**

Una de las principales fuentes de datos raster a nivel mundial es *EARTH DATA* de la NASA y puedes descargar los raster en [https://disc.gsfc.nasa.gov](https://disc.gsfc.nasa.gov).

- **Librería getSpatialData**

La librería getSpatialData provee raster de los proyetos Sentinel-1, Sentinel-2, Sentinel-3, Sentinel-5P, Landsat 8 OLI, Landsat ETM, Landsat TM, Landsat MSS, MODIS (Terra & Aqua) y SRTM DEMs de la NASA. Puede consultar el repositorio de la libreria en [https://github.com/16EAGLE/getSpatialData](https://github.com/16EAGLE/getSpatialData).

---
class: inverse, center, middle
name: read
# [2. Importar datos raster](#contenido)

---
# 2.1. Cargemos los datos (...)

```r
# llamar las librerias
require(raster)

# Leer raster
deforestacion = raster("data/input/siac/magdalena_deforestacion_1990_2000.tif")
```

```
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum Unknown based on GRS80 ellipsoid in Proj4
## definition
```

```r
# atributos del raster
deforestacion
```

```
## class      : RasterLayer 
## dimensions : 2505, 2704, 6773520  (nrow, ncol, ncell)
## resolution : 0.0002753217, 0.0002713948  (x, y)
## extent     : -74.28793, -73.54346, 10.42546, 11.1053  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +no_defs 
## source     : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/SIAC/magdalena_deforestacion_1990_2000.tif 
## names      : magdalena_deforestacion_1990_2000 
## values     : 1, 5  (min, max)
```

---
# 2.2. Veamos el raster

```r
# plot basico del raster
plot(deforestacion)
```
<div align="center">
<img src="pics/Rplot.png" height=450>
</div>

---
# 2.3. Geometría del raster

```r
# nombres de las bandas
names(deforestacion) 
```

```
## [1] "magdalena_deforestacion_1990_2000"
```

```r
names(deforestacion) = "cobertura_vegetal"
names(deforestacion) 
```

```
## [1] "cobertura_vegetal"
```

```r
# Extension
st_bbox(deforestacion)
```

```
##      xmin      ymin      xmax      ymax 
## -74.28793  10.42546 -73.54346  11.10530
```

```r
# Proyeccion
deforestacion@crs
```

```
## CRS arguments: +proj=longlat +ellps=GRS80 +no_defs
```

```r
crs(deforestacion)
```

```
## CRS arguments: +proj=longlat +ellps=GRS80 +no_defs
```

---
# 2.4. Valores del raster (...)

```r
# Valores de los pixeles
minValue(deforestacion$cobertura_vegetal)  
```

```
## [1] 1
```

```r
maxValue(deforestacion$cobertura_vegetal)  
```

```
## [1] 5
```

```r
values(deforestacion) %>% head()
```

```
## [1] NA NA NA NA NA NA
```

```r
# Cambiar los values de los pixeles
values(deforestacion)[1] <- 0
values(deforestacion)[is.na(values(deforestacion))==T] <- 0
values(deforestacion) %>% head()
```

```
## [1] 0 0 0 0 0 0
```

---
# 2.4. Valores del raster (...)

```r
# Descriptivas
values(deforestacion) %>% table()
```

```
## .
##       0       1       2       3       4       5 
##   76469 2432531  200575   88669   26717 3948559
```

```r
values(deforestacion) %>% summary()
```

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   1.000   5.000   3.388   5.000   5.000
```

```r
value_raster = values(deforestacion) 
str(value_raster)
```

```
##  num [1:6773520] 0 0 0 0 0 0 0 0 0 0 ...
```

---
class: inverse, center, middle
name: ope
# [3. Operaciones geométricas con raster](#contenido)

---
# 3.1. Importar raster de luces

```r
# leer raster
luces = raster(x = 'data/input/NOAA/colombia_202004.tif')
luces
```

```
## class      : RasterLayer 
## dimensions : 2989, 2918, 8721902  (nrow, ncol, ncell)
## resolution : 0.004166667, 0.004166667  (x, y)
## extent     : -79.00625, -66.84792, 0.002082733, 12.45625  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NOAA/colombia_202004.tif 
## names      : colombia_202004 
## values     : -0.43, 1389.84  (min, max)
```

```r
# leer polygono
cataca = st_read(dsn= 'data/input/mgn/MGN_Municipio.shp') %>% subset(MPIO_CCDGO %in% c('053'))
```

```
## Reading layer `MGN_Municipio' from data source `/Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/mgn/MGN_Municipio.shp' using driver `ESRI Shapefile'
## Simple feature collection with 30 features and 12 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.9471 ymin: 8.912502 xmax: -73.54336 ymax: 11.34912
## geographic CRS: WGS 84
```

```r
cataca
```

```
## Simple feature collection with 1 feature and 12 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.25773 ymin: 10.42533 xmax: -73.54336 ymax: 10.8791
## geographic CRS: WGS 84
##   OBJECTID DPTO_DPTO_ MPIO_CCDGO MPIO_CNMBR
## 3      659         47        053  ARACATACA
##                             MPIO_CRSLC MPIO_NANO_ MPIO_NAREA MPIO_CSMBL
## 3 Ordenanza 47 del 28 de Abril de 1915       1993   1746.793          4
##   MPIO_CCNCT MPIO_NANO SHAPE_Leng SHAPE_Area                       geometry
## 3      47053      2012   2.650507  0.1443276 POLYGON ((-73.72812 10.8626...
```

---
# 3.2. Cliping un raster(...)

```r
l_cataca = crop(luces,cataca) # Qué debo revisar primero?
ggplot() + geom_tile(data = as.data.frame(l_cataca, xy=TRUE), aes(y=y,x=x,fill=colombia_202004)) + 
scale_fill_distiller(palette='Spectral',na.value = 'gray') + 
geom_sf(data =cataca,color = 'black',fill=NA) + theme_bw()
```

---
# 3.2. Cliping un raster(...)

```r
l_cataca = crop(luces,cataca) %>% mask(cataca)
ggplot() + geom_tile(data = as.data.frame(l_cataca, xy=TRUE), aes(y=y,x=x,fill=colombia_202004)) + 
scale_fill_distiller(palette='Spectral',na.value = 'gray') + 
geom_sf(data =cataca,color = 'black',fill=NA) + theme_bw()
```

---
# 3.3. Extraer los valores del raster

```r
values(l_cataca) %>% .[is.na(.)==F] %>% head()
```

```
## [1] 0.00 0.00 0.00 0.00 0.21 0.00
```

```r
data = values(l_cataca) %>% .[is.na(.)==F]
summary(data)
```

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.3600  0.2600  0.3300  0.3862  0.4000 19.7600
```

---
# 3.4. Raster a datos vectoriales (...)

```r
# Raster a puntos
point = rasterToPoints(l_cataca,spatial = T) %>% st_as_sf()
class(point)
```

```
## [1] "sf"         "data.frame"
```

```r
ggplot() + geom_sf(data = point , aes(color=colombia_202004),size=0.1)
```

---
# 3.4. Raster a datos vectoriales (...)

```r
# Raster a polygonos
polygon = rasterToPolygons(l_cataca) %>% st_as_sf()
ggplot() + geom_sf(data = polygon , aes(fill=colombia_202004),size=0) +
scale_fill_distiller(palette='Spectral',na.value = 'gray')
```

---
class: inverse, center, middle
name: layers
# [4. Raster de varios layers](#contenido)

---
# 4.1. Cómo se ve un raster RGB?
<div align="center">
<img src="pics/rgb_raster.png" height=500>
</div>
Tomado de: [https://www.neonscience.org](https://www.neonscience.org)

---
# 4.2. Importar raster

```r
GDALinfo("data/input/neon/HARV_RGB_Ortho.tif")
```

```
## rows        2317 
## columns     3073 
## bands       3 
## lower left origin.x        731998.5 
## lower left origin.y        4712956 
## res.x       0.25 
## res.y       0.25 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## file        data/input/neon/HARV_RGB_Ortho.tif 
## apparent band summary:
##    GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64           TRUE   -1.7e+308          1       3073
## 2 Float64           TRUE   -1.7e+308          1       3073
## 3 Float64           TRUE   -1.7e+308          1       3073
## apparent band statistics:
##   Bmin Bmax Bmean Bsd
## 1    0  255   NaN NaN
## 2    0  255   NaN NaN
## 3    0  255   NaN NaN
## Metadata:
## AREA_OR_POINT=Area
```

---
# 4.3. Veamos que tenemos

```r
banda_r <- raster(x = "data/input/neon/HARV_RGB_Ortho.tif") # R carga por 'default' la banda 1, es decir la roja
paleta_col <- gray.colors(n = 100, start = 0.0,end = 1.0,alpha = NULL) 
plot(banda_r, col=paleta_col, axes=FALSE, main="Imagen RGB - Banda 1 (roja)") 
```

---
# 4.4. Values de la banda
En un RGB podemos tener 255*255*255 posibles combinaciones, es decir 16.581.375 colores.

```r
minValue(banda_r)
```

```
## [1] 0
```

```r
maxValue(banda_r)
```

```
## [1] 255
```

```r
values(banda_r) %>% summary()
```

```
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0    81.0   108.0   102.2   128.0   255.0
```

---
# 4.5. cargar las 3 bandas

```r
RGB_stack = stack("data/input/neon/HARV_RGB_Ortho.tif")
plotRGB(RGB_stack, r = 1, g = 2, b = 3)
```

---
# 4.6. Inspeccionemos el objeto (...)

```r
# tipo de objeto
class(RGB_stack)
```

```
## [1] "RasterStack"
## attr(,"package")
## [1] "raster"
```

```r
dim(RGB_stack)
```

```
## [1] 2317 3073    3
```

```r
# nombres de las bandas
names(RGB_stack)
```

```
## [1] "HARV_RGB_Ortho.1" "HARV_RGB_Ortho.2" "HARV_RGB_Ortho.3"
```

```r
names(RGB_stack) = c('red','green','blue') 
names(RGB_stack)
```

```
## [1] "red"   "green" "blue"
```
---
# 4.6. Inspeccionemos el objeto (...)

```r
RGB_stack@layers # ver las bandas
```

```
## [[1]]
## class      : RasterLayer 
## band       : 1  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NEON/HARV_RGB_Ortho.tif 
## names      : red 
## values     : 0, 255  (min, max)
## 
## 
## [[2]]
## class      : RasterLayer 
## band       : 2  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NEON/HARV_RGB_Ortho.tif 
## names      : green 
## values     : 0, 255  (min, max)
## 
## 
## [[3]]
## class      : RasterLayer 
## band       : 3  (of  3  bands)
## dimensions : 2317, 3073, 7120141  (nrow, ncol, ncell)
## resolution : 0.25, 0.25  (x, y)
## extent     : 731998.5, 732766.8, 4712956, 4713536  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : /Users/eduard-martinez/Dropbox/teaching/Taller de R/GitHub/clases/clase-13/data/input/NEON/HARV_RGB_Ortho.tif 
## names      : blue 
## values     : 0, 255  (min, max)
```

---
### 4.7. Extraer atributos

```r
point_RGB = rasterToPoints(RGB_stack,spatial = T) %>% st_as_sf()
point_RGB
```

```
## Simple feature collection with 7120141 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 731998.6 ymin: 4712956 xmax: 732766.6 ymax: 4713535
## CRS:            +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## First 10 features:
##    red green blue                 geometry
## 1    0     1    0 POINT (731998.6 4713535)
## 2    2     0   10 POINT (731998.9 4713535)
## 3    6     9    1 POINT (731999.1 4713535)
## 4    0     0    0 POINT (731999.4 4713535)
## 5   16     5   17 POINT (731999.6 4713535)
## 6    0     0    0 POINT (731999.9 4713535)
## 7    0     4    4 POINT (732000.1 4713535)
## 8    6     2    0 POINT (732000.4 4713535)
## 9    1     1    0 POINT (732000.6 4713535)
## 10   5     0    7 POINT (732000.9 4713535)
```

---
class: inverse, center, middle
name: view

# [Hoy vimos...](#contenido1)
<html><div style='float:left'></div><hr color='#FF007F' size=1px width=796px></html>

---
#  Hoy vimos...

- ☑ [Introducción a datos raster](#intro)

- ☑ [Importar datos raster](#read)

- ☑ [Operaciones geométricas con raster](#ope)

- ☑ [Raster de varios layers](#layers)